# install.packages("rgbif")
data <- rgbif::occ_data(scientificName = "Jubaea chilensis",
limit=10000,
country = "CL",
basisOfRecord = "HUMAN_OBSERVATION")$dataAnálisis Espacial
Minicurso SOCHE
1 Análisis de patrones de puntos
1.1 Cargar datos
Trabajaremos con un set de datos que descargaremos con el paquete {rgbif}.
GBIF —Infraestructura Mundial de Información en Biodiversidad— es una organización internacional y una red de datos financiada por gobiernos de todo el mundo, destinada a proporcionar a cualquier persona, en cualquier lugar, acceso abierto y gratuito a datos sobre cualquier tipo de forma de vida que hay en la Tierra.
ESto está en formato Darwin Core:
El set de datos también se encuentra disponible en el github.
1.2 Convertir a objeto espacial
Al descargar este set de datos, estará disponible como un data.frame por lo que es necesario transformarlo a un objeto espacial para hacer algunas operaciones. Esto lo haremos con el paquete {sf}:
# install.packages("sf")
# install.packages("tidyverse")
library(sf)
library(tidyverse)
data_sp <- data %>%
filter(!is.na(decimalLongitude), !is.na(decimalLatitude)) %>%
st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)1.3 Visualizar
Podemos visualizarlo de manera rápida con el paquete {mapview}:
# install.packages("mapview")
library(mapview)
mapview(data_sp)Ahora podemos comenzar con nuestro análisis de patrones de puntos:
1.4 Corroborar sistema de proyección
Este tipo de análisis solo puede ser hecho con sistemas de coordinadas proyectados, entonces es necesario cambiar nuestro sistema de coordenadas:
data_sp_utm <- data_sp %>%
st_transform(crs = 32719)
data_sp_utm$geometryGeometry set for 360 features
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 221553.9 ymin: 6075129 xmax: 356415 ymax: 6692283
Projected CRS: WGS 84 / UTM zone 19S
First 5 geometries:POINT (259482.3 6337580)POINT (259505 6337556)POINT (259517.3 6337518)POINT (341457.1 6236730)POINT (258941.6 6337451)st_crs(data_sp_utm)Coordinate Reference System:
User input: EPSG:32719
wkt:
PROJCRS["WGS 84 / UTM zone 19S",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 19S",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-69,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",10000000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["Between 72°W and 66°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Bolivia. Brazil. Chile. Colombia. Peru."],
BBOX[-80,-72,0,-66]],
ID["EPSG",32719]]st_crs(data_sp)Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]1.5 Establecer ventana de análisis
Como primera aproximación utilizaré la extensión de los puntos:
win <- data_sp_utm %>%
st_bbox() %>%
st_as_sfc()
mapview(list(data_sp_utm, win))1.6 Crear patrón de puntos
Ahora utilizaremos el paquete {spatstat} para realizar nuestro análisis. Para crear un patrón de puntos es necesario tener las coordenadas de los puntos y la ventana de análisis:
# install.packages("spatstat")
library(spatstat)
x <- st_coordinates(data_sp_utm)[,1]
y <- st_coordinates(data_sp_utm)[,2]
box <- as.owin(win)
pp1 <- ppp(x = x, y = y, window = box)
plot(pp1)
1.7 Realizar kernel de densidad
Para calcular la intensidad de uso utilizaremos la función density:
kernel <- density.ppp(pp1)
plot(kernel)
library(stars)Warning: package 'stars' was built under R version 4.2.1Loading required package: abindkernel_stars <- kernel %>% st_as_stars() %>% st_set_crs(32719)
mapview(kernel_stars)Podemos estandarizar el kernel haciendo un poco de álgebra:
kernel_stars$valor_01 <- (kernel_stars$v - min(kernel_stars$v))/(max(kernel_stars$v)- min(kernel_stars$v))
mapview(kernel_stars[2])1.8 K de ripley
La función K de Ripley, K(r) es un método basado en la distancia que mide la aglomeración de un patrón de puntos espacial contando el número medio de vecinos que presenta cada punto dentro de un círculo de radio (r) determinado, en un determinado espacio, la función K compara el valor observado a una cierta distancia con el valor esperado a esa misma distancia; dado un proceso de Poisson homogéneo, también conocido como complete spatial randomness (CSR), es decir, que todos los puntos tienen la misma probabilidad de ocurrir en cualquier parte del área de estudio Fuente
Hacemos el análisis de K de Ripley con la función Kest
k_ripley <- Kest(pp1)
plot(k_ripley)
Para más información sobre la interpretación, VER